% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% This code produces  Table 3 Panel D
clear;
tic
cd  ..\..\Simulations_VC\MS_produits_KR_CalvoPlus

load  ..\..\Simulations_VC\MS_produits_KR_CalvoPlus\stat4_alt
param=stat4_alt;
param(143,:)=[0.857199875280086 0.266695001309508 0.0297130383012599 0.694600000000000 0.743983380443812 0.529997861050279 0.0254200000000000 0.0508400000000000 2.77502486296478 0.0139400000000000 0.0254200000000000 0.0393600000000001 1 0.200021250982131 0 0.000693548960074729]


load  ..\..\Estimation_param\MS_produits_KR_CalvoPlus\actual_moments_k
weight=actual_moments_k(:,12);
sum(weight);
secteur=actual_moments_k(:,1);

p0=param(:,1);
%muc=param(:,2);
sig_ea=param(:,3);
p_a=param(:,4);
%sig_a=sig_ea./sqrt(1-rho_a.^2);
share_calvo=param(:,end-3);
muc=param(:,end-2);
av_muc=param(:,end-1);

sig_tilda=sqrt((2*p_a.*sig_ea.^2)/(1+0.6946))
mu_tilda=muc.*(1-share_calvo);


stat=[p0 muc mu_tilda sig_ea p_a  sig_tilda share_calvo];
%mean(stat)
%median(stat)
size_param=size(stat,2)

stat_q1=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q1(i) = mean(sortx([j j+1]));
            else
                
                stat_q1(i) = sortx(j+1);
            end
     stat_q1(end)= csumw(end);
     stat_q1(end-1)= v;
    end  
 end


stat_med=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_med(i) = mean(sortx([j j+1]));
            else
                
                stat_med(i) = sortx(j+1);
            end
     stat_med(end)= csumw(end);
     stat_med(end-1)= v;
    end  
 end
%"stat_med"
%stat_med


stat_q3=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = 3*sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q3(i) = mean(sortx([j j+1]));
            else
                
                stat_q3(i) = sortx(j+1);
            end
     stat_q3(end)= csumw(end);
     stat_q3(end-1)= v;
    end  
 end



meanw_data_all=weight'*stat/sum(weight);
for i=1:size(stat,2)
stdw_data(i)=sqrt(weight'*((stat(:,i)-meanw_data_all(i)).^2)/sum(weight));
end
%fprintf("Stats table Estimates - PANEL D")
[meanw_data_all(1:size_param);
stat_med(1:size_param);
stat_q3(1:size_param)-stat_q1(1:size_param);
stdw_data(1:size_param)]
%'interq / med'
%(stat_q3-stat_q1)./meanw_data_all
%(stat_q3-stat_q1)./stat_med


%IDEM Excluding energy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stat=stat((stat(:,1)<=0.3),:);
weight=weight((stat(:,1)<=0.3),:);

stat_q1=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q1(i) = mean(sortx([j j+1]));
            else
                
                stat_q1(i) = sortx(j+1);
            end
     stat_q1(end)= csumw(end);
     stat_q1(end-1)= v;
    end  
 end


stat_med=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_med(i) = mean(sortx([j j+1]));
            else
                
                stat_med(i) = sortx(j+1);
            end
     stat_med(end)= csumw(end);
     stat_med(end-1)= v;
    end  
 end

stat_q3=ones(1,size_param+2);
stat_sect=stat;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = 3*sum(sortw)/4;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_q3(i) = mean(sortx([j j+1]));
            else
                
                stat_q3(i) = sortx(j+1);
            end
     stat_q3(end)= csumw(end);
     stat_q3(end-1)= v;
    end  
 end


meanw_data_all=weight'*stat/sum(weight);
for i=1:size(stat,2)
stdw_data(i)=sqrt(weight'*((stat(:,i)-meanw_data_all(i)).^2)/sum(weight));
end

%fprintf("Idem excluding energy")

[meanw_data_all(1:7);
stat_med(1:7);
stat_q3(1:7)-stat_q1(1:7);
stdw_data(1:7)]
'interq / med'
(stat_q3-stat_q1)./stat_med


